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We propose a hybrid deterministic and stochastic approach to achieve extended time scales in 
atomistic simulations that combines the strengths of Molecular Dynamics (MD) and Monte Carlo 
(MC) simulations in an easy-to-implement way The method exploits the rare event nature of the 
dynamics similar to most current accelerated MD approaches but goes beyond them by providing, 
without any further computational overhead, (a) rapid thermalization between infrequent events, 
thereby minimizing spurious correlations and (b) control over accuracy of time scale correction, 
while still providing similar or higher boosts in computational efficiency We present two applica- 
tions of the method: (a) vacancy mediated diffusion in Fe yields correct diffusivities over a wide 
range of temperatures and (b) source controlled plasticity and deformation behavior in Au nanopil- 
lars at realistic strain rates (10^/sec and lower) with excellent agreement with previous theoretical 
predictions and in situ high-resolution transmission electron microscopy (HRTEM) observations. 
The method gives several orders of magnitude improvements in computational efficiency relative to 
standard MD and good scalability with size of system. 



With vast improvements in the quality of available in- 
teratomic force-fields and computer power, the classical 
MD simulation has seen a dramatic increase in its use 
across a variety of fields over the past few decades 
One of the features that makes MD so appealing is its 
ability to actually follow the dynamical evolution of the 
system, thus giving insight into the microscopic behavior 
of the material. However, this is where the major limita- 
tion of MD comes into light too: most of the interesting 
dynamics occurs as the system moves from one energy 
basin to another through infrequent rare events, while the 
system remains stuck in some energy basin for extended 
periods of time. This non-ergodicity, coupled with the 
small time steps (on the order of femtoseconds) needed 
for total energy staying conserved, severely restricts the 
timescales accessible in MD simulations and also leads 
to limited phase space exploration. There have been 
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FIG. 1: (Color online) Schematic and flowchart of the algo- 
rithm. Shell Sw of constant potential energy and energy well 
W as described in text are shown here. See Eq.® for defini- 
tion of ?w- 

many attempts at addressing this time-scale problem in 
MD - examples include metadynamics^, the activation 



relaxation technique^, parallel replica dynamics, temper- 
ature accelerated dynamics and hyperdynamic^*^. There 
are several excellent reviews such as Ref. ^ available 
on the subject. The hyperdynamics method^ offers an 
elegant and practical way to increase the rate of infre- 
quent events. It consists of adding a potential energy bias 
that makes the potential wells, in which the system nor- 
mally remains trapped for extended periods, less deep. 
A time-scale correction is also evaluated in terms of the 
bias potential. The hyperdynamics method, especially 
with the advent of a variety of easy to implement biasing 
forms^, has seen several compelling applications over the 
past years^— . Our approach in this letter builds upon 
the crucial insights of Voter and co-workers while seeking 
improvements along two important dimensions. First, it 
bypasses a fundamental trade-off present in hyperdynam- 
ics: a shallower potential well provides faster dynamics 
but, at the same time, reduces the ability of the modified 
potential to properly thermalize the system in between 
the infrequent events, resulting in artificial correlation 
between these events. Second, our method provides bet- 
ter independent control over the accuracy of time scale 
correction, while hyperdynamics time scale estimates can 
remain noisy up to long simulation times, especially for 
large system sizes (see Ref. 14 for a discussion on this). 

Let the state of the system be characterized by position 
X and velocity each being a 3 A/"- dimensional vector for 
a system of A^ atoms. When the potential energy V{x) 
of the system is above a threshold Vq , the system evolves 
via constant-energy (or constant-temperature) MD ac- 
cording to its true Hamiltonian (Fig. [1]). This high 
energy region of the phase space contains the interest- 
ing but infrequently occuring events. The method is for- 
mally correct for any choice of l^; a higher choice of Vq 
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merely limits our ability to monitor the detailed dynam- 
ics of some events. When the system's potential energy 
falls below Vq^ two MC simulations are initiated (denoted 
a and h). Simulation a runs MC with a perfectly uni- 
form potential inside the potential well W consisting the 
points X where the true potential energy V{x) lies be- 
low Vo (i.e. all moves are accepted as long as they do 
not go outside the well). Simulation a is run until the 
system is well thermalized and has lost memory of how 
it entered the well (this takes a few MC passes, an in- 
significant amount of wall clock time). MD then resumes 
with positions drawn from the last MC state that visited 
the boundary of the potential well. The vector v of the 
velocities of all atoms for restarting MD is drawn from 
a Maxwell-Boltzmann distribution corresponding to the 
temperature T of interest, conditional on v ■ \/V{x) > 
(i.e. we only consider velocities in the half-space point- 
ing outwards of the well). MC simulation a is first of 
the crucial differences between our approach and hyper- 
dynamics: it ensures proper thermalization of the system 
between rare events even when using a completely fiat po- 
tential in the well. Even though it is done with a uniform 
potential, it does not lead to the molecular structure be- 
ing completely lost since we rule out all moves that lead 
to energy higher than Vq. 

In parallel to simulation a, another MC simulation b 
is launched to estimate the mean time the system should 
have spent in the well W. Akin to simulation a, b also 
rejects all moves that land outside the well W. The mean 
time spent in W is given by the reciprocal of the flux 
exiting^^ the well W: 
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where the average (• • • ) is taken over x drawn from 
the well W with a probability density proportional to 
^-v{x)/{kBT) y^i^QYQ j^^ Boltzmann's constant and the 
following definitions hold: 1{A) equals 1 if the event A 
is true and otherwise, S^^ is a shell of width w at the 
boundary of the well W, which can be defined in the 
limit of small w as 



S^ = {x: \V{x) - Vol < w\\/V{x)\/2} 



(2) 



and V denotes the mean projection of a Maxwell- 
Boltzmann-distributed velocity along the unit vector u 
parallel to VF(x), conditional on v ■ u > 0. The latter is 
given by 
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where rrii is the mass of atom i and \ui \ denotes length of 
the 3 dimensional subvector of u associated with atom i. 
Note that the Eq.(|3]) reduces to the familiar expression^^ 
V = Y^/c^T/27rm when all atoms have the same mass, 
in which case v factors out of the average in ([1]). Since 
Eq.([T]) involves an average, it can be approximated by 



MC simulations. However, the most straightforward im- 
plementation of this approach would be very inefficient 
because x would rarely visits the boundary of the 
well. The efficiency can be considerably improved by us- 
ing a biased potential V*{x) which is the same as the real 
potential in the high energy regions (i.e. regions outside 
well W with V{x) > Vb), but lifted up in the deep energy 
basins. With this Eq. ([T]) becomes 
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where the averages (•••)* are taken over x drawn from 
the well W with a probability density proportional to 
^-v*ix)/ikBT) p ^ l/ksT. MC simulation b is the 
second main difference with hyperdynamics: it provides 
separate control over the accuracy of the speed up fac- 
tor since the the length of the MC simulation b can be 
adjusted independently of the length of the whole simu- 
lation. 

The form of biasing we use is a well established and 
easy to implement biasing potential used in several imple- 
mentations of Voter's hyperdynamics method, proposed 
by Hamelberg et 
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The times tw obtained via MC simulations b can be 
directly added to the physical time spent doing MD sim- 
ulations to yield the overall physical time of the simu- 
lation. However, refinements of the method can yield 
further improvement in efficiency. Instead of computing 
for each well W, one may keep a running average 
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of the time spend in the rib wells sampled via MC simula- 
tion b {ua^ the number of wells actually visited, may well 
far exceed n^). Once this average is converged, there is 
no need to initiate MC simulation b anymore. The over- 
all time spent in all the wells will simply be tw "^ria/nh. 
Note that there is no need to keep separate averages for 
different types of wells, which would have been difficult 
to implement. Although MC simulations a still need to 
be performed for all wells, these converge much more 
rapidly. Other efficiency improvements can be obtained 
by not performing fully converged MC simulations b and 
exploiting the fact that errors will average out over wells 
in Eq.(j6]). Note that this scheme must be used while 
ensuring that the biasing potential is sufficiently strong 
so that most of the random errors in Eq.Q are concen- 
trated in the numerator, to avoid a systematic bias due 
to nonlinear ity of the ratio. We would like to point out 
that only the parameter w is additional to those in any 
typical hypderdynamics scheme (Hamelberg et al.^^s in 
this case), the choice of which does not effect the re- 
sult since we extrapolate tw to the limit of small w^. 
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FIG. 2: (Color online) Loge (diffusivity) versus inverse tem- 
perature, (a) Straight lines denote Mendelev et al^SF^ MD 
calculations. These are valid only until 700 K. (b) Aster- 
isks denote diffusivity measurements per our approach, (c) 
The dashed line shows experimental measurements^^ valid be- 
tween 1000 and 1200 K. 



Our approach compares favorably with hyperdynamics^ 
where one does not have control over the accuracy of the 
accelerated time (hyperdynamics relies on this error can- 
celling out over time but does not provide an estimate of 
how much it is^ii^), and one is obliged to keep perform- 
ing dynamics with the biased potential at all stages of 
the calculation. Thus our method offers boosts as high 
as those that one could get from setting a = in Eq.(|5]) 
(akin to the flooding scheme^^^^^), but still avoiding the 
slow convergence in time and problems with discontinu- 
ous forces that one encounters in implementing flooding 
based hyperdynamics. In addition we avoid errors from 
sampling the system in the state when it is not thermal- 
ized between two rare events - once MD is relaunched in 
our scheme, the system is well thermalized by virtue of 
simulation a. To minimize the wall-clock time needed for 
calculation of time in Eq.(j4]) via simulation 6, we use an 
optimal extent of biasing as suggested in Ref. 9. This 
involves setting a o:^ Vq — Vmin which allows the biased 
potential to capture the shape of the potential wells^. a 
smaller than this would improve sampling of the numera- 
tor in Eq.dH) but detoriate that of the denominator. We 
picked two problems to demonstrate that our method 
yields correct dynamics: (a) vacancy diffusion in BCC 
Fe at room temperature, and (b) deformation behavior 
in Au nanopillars at realistic strain rates. 

Lattice diffusion at low temperatures is beyond the 
time scales one can access in current MD simulations, 
with most investigations^^ only beyond 700 K. The sys- 
tem we consider is 249 Fe atoms (5x5x5 BCC supercell 
with 1 vacancy) interacting through the Embedded Atom 
Method (EAM) potential^^. For the MD part here and in 
deformation behavior problem, we performed NVT sim- 
ulations using time step of 2xl0~^^ sec and a Langevin 



FIG. 3: (Color online) Simulation cell for stress-strain calcula- 
tions, (a) Prior to application of any strain, (b) after yielding 
(strain = 12%) with strain rate = 5xl0^/sec. In (b) the lead- 
ing partial has nucleated on {111} slip plane leaving behind 
the 2-layer thick HCP region denoting an intrinsic stacking 
fault. Failure is thus through slip and not twinning, in agree- 
ment with Ref. 25. Atoms are identified as per bond order 
parameter Qe^^'^. Perfect HCP atoms have been removed 
for clarity. 



thermostat with coupling constant 1x10 sec ^. The 
biasing parameter a was 50 eV. The Vq values we used 
at 500 and 300 K were -975.5 eV and -984 eV resp. (4 
and 2.5 eV more than the mean energy at 500 and 300K 
resp.). We took the equilbrium concentration of defects^^* 
to convert our effective diffusivity into equilibrium dif- 
fusivity. In Fig. [2] we plot the equilibrium diffusivity 
as obtained from (a) MD simulations ^^,(b) our proposed 
approach, and (c) experimental measurement^ that in- 
clude contributions from interstitial migrations also and 
hence are only slightly higher than both ours and MD 
values. We obtain around 5 orders of magnitude boost, 
with similar speed up factors for system sizes up to 30000 
atoms. 

For our second problem (see Fig. [3]), we looked at 
deformation behavior of Au nanopillars. With advent of 
excellent in situ TEM and HRTEM tools, there are many 
elegant experiments on sub-lO-nm sized crystals^^^^^^^^. 
Deformation in such small sizes is controlled by dis- 
location nucleation, and has been phenomenologically 
predicted^^ and experimentally found— i^^i^S to have 
small activation volumes and strong strain-rate sensitiv- 
ity. However there is no direct MD based confirmation of 
this strong strain-rate sensitivity due to inability of MD 
to reach strain rates lower than 10^/sec. 

Using our method we were able to reach 10^/sec strain- 
rate regime with only around 48 hours of computer 
time. We could also obtain several correct qualitative 
and quantitative aspects of the deformation dynamics, 
without assuming anything about the nature of deforma- 
tion. The system we consider is 2016 Au atoms (cylin- 
der with height 7.4 nm and diameter 2.5 nm) interact- 
ing through EAM potential^^. The biasing parameter 
a was 1500 eV while the starting Vb value used was 
-7266 eV(8 eV more than the mean energy at 300K), 
adjusted every 1000 MD steps to take into account the 
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FIG. 4: (Color online) (a) Stress-strain plots at 3 different 
strain rates: 2.5x10^ /sec (open circles), 5x10^ /sec (asterisks), 
5x10^ /sec (pluses). The initial stress corresponding to zero- 
strain is non-zero due to surface effects^, (b) loge (stress) at 
11% strain (relative to surface stress at zero strain) versus 
loge (strain rate). 

pressure- volume work contribution from the stress. The 
cylinder was initially carved out from perfect FCC lat- 
tice (Fig. [3l^a)). Periodic boundary conditions were im- 
posed only along the cylinder axis z which is the same 
as the compression axis (001). The cylinder was first 
equilibrated for 500 ps before beginning the compression 
carried out by uniformly re-scaling the z-coordinates of 
all atoms. The atomic virial stress^^ was used to obtain 
the Cauchy stress. 4 different strain rates e were consid- 
ered: 5xl0^/sec, 2.5xl0^/sec, SxlO^sec and SxlO'^/sec. 
We present the resulting stress(cr)-strain(e) plots in Fig. 
[11(a). Several conclusions can be drawn from Figs. [3] and 
[4] that prove our algorithm capable of predicting correct 
dynamics in complicated systems. The yielding occurs 
around 10% strain, and is through slip and not twinning 



or elastic instabilities: a leading partial nucleates on a 
{111} slip plane at lower stresses than a trailing partial. 
This can be seen in Fig. [31(b) where the leading partial 
nucleated from the surface and left behind a 2-layer thick 
HCP region which again changes back to FCC after the 
trailing partial also nucleates at higher stresses and re- 
combines with the leading partial. Fig. [Sfb) is identical 
to HRTEM images for (001) loading of Au nanowires^^. 
The strain rate sensitivity m in the relation a = <Joe^ 
(derived by looking at stress at 11% strain) is around 
0.14±0.07 (see Fig. [^b), while Ref. 20 reports it to be 
around 0.11 for 75 nm diameter pillars. The activation 
volume ^ for the deformation as calculated through^^« 
= y/?>kBTd{lge) / d(j\^ around Ih^ (6=burgers vector) in 
excellent agreement with experiments observations^flaSi. 

To summarize, we have proposed an approach that 
combines the strengths of MC and MD, thus offering 
boosts of several orders of magnitudes with good system 
size scaling. We have applied the method to study lattice 
diffusion in BCC Fe at low temperatures and deforma- 
tion of Au nanopillars at low strain rates and found it to 
work really well in both cases, predicting correct dynam- 
ics and exhibiting good scaling with increase in system 
size from 249 to 2016 atoms. We thus expect the method 
to be useful in a variety of situations. 
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